1 Goal

This notebook can be called with knit_notebook8_byESM.R; it can be called in a loop over ESMs and target experiments, and will knit an analysis html for each.

These values must be defined globally in knit_notebook8_byESM.R or uncommented and defined here to run this notebook directly:

# # what we will emulate:
# esm_name <- "MPI-ESM1-2-LR"
# esm_experiment <- 'ssp245'
# 
# # How many times do we want to draw full generated ensembles:
# Ndraws <- 1

print(esm_name)

## [1] "CanESM5"

print(esm_experiment)

## [1] "ssp245"

We will only explore tol=0.1 degC for matching in this notebook (per CT, from the SD of the smoothed Tgav time series. Should check for each ESM, for now just use 0.1 value from CanESM).

2 Set Up

# Then load the functions we will want to use, these are currently written in R and will be translated into python. 
source("functions/nearest_neighbor_matching.R") # the function match_neighborhood is defined here
source("functions/permuting_matches.R")  # the function permute_stitching_recipe is defined here
source("functions/stitching_functions.R")  # the function stitch_global_mean is definend here 
source("functions/gridded_recipe.R") # wrappers for the work flow from matching to output csv needed by python for 
                                     # layering in gridded data to output as netcdf

set.seed(1)

Load the inputs that will be used. In addition to SSP534-over having enough data issues to fully exclude, at least some of the realizations for SSP434 and SSP460 just don’t have future data at all for CanESM5. And for other models, some realizations have odd windows (Like a 2012-2014 window) and those realizations are excluded as well. We need to go back to pangeo and fix this for sure, but for now, I’ll work with the archive I have. Because we don’t want to lose too much of the archive now, I will only remove those realizations and not a full experiment for a single bad data point.

# Time series of the raw global mean temperature anomaly, this is the data that is 
# going to be stitched together. Historical data is pasted into every scenario. No
# smoothing has been done.
tgav_data <- read.csv("inputs/main_raw_pasted_tgav_anomaly_all_pangeo_list_models.csv", 
                      stringsAsFactors = FALSE) %>% dplyr::select(-X)


# A chunked smoothed tgav anomaly archive of data. This was tgav_data but several 
# steps were taken in python to transform it into the "chunked" data and save it 
# so that we do not have to repeate this process so many times. 
read.csv(paste0('inputs/', esm_name, '_archive_data.csv'), stringsAsFactors = FALSE) %>%
  # Exclude ssp534-over:
  dplyr::filter(experiment != 'ssp534-over') %>%
  # eliminate any individual realizationXexperiment runs that don't have all the years:
  group_by(experiment, variable,   model  , ensemble) %>%
  mutate(nWindows = n()) %>%
  ungroup %>%
  filter(nWindows == 28) %>%
  select(-nWindows) ->
  archive_data


# The target data - all ssp245 realizations (smoothed because it's the target for
# matching)
archive_data %>%
  filter(experiment == esm_experiment) ->
  target_data


# The corresponding unsmoothed original data for validation
tgav_data %>%  
  filter(model == unique(target_data$model) & experiment == unique(target_data$experiment))  %>%
  select(-file, -timestep, -grid_type) -> 
  original_data

We also need to prepare the file of path metadata for pangeo. This is a csv file that contains the pangeo paths to different variables of interest. By selecting the paths for tas, psl, and pr, keeping only those experimentXensemble combinations that have defined paths to netcdfs for all three variables, and then using the resulting experiment X ensemble combinations as an additional filter on the archives we use for matching, we make sure that we can then stitch together gridded data for all 3 of these variables no matter what our generated ensembles are. Adding or removing variables of interest is relatively trivial. For now, the path metadata file has information for tas, psl, pr, tasmin, and tasmax.

# The main pangeo file of path metadata so that we can filter the archives to include only
# the experiment*ensemble members with netcdfs for all variables of interest
read.csv("inputs/pangeo_path_metadata_tas_psl_pr_tmax_tmin.csv", 
         stringsAsFactors = FALSE,
         na.strings = c("", " ", "NA")) %>%  
  rename(model = source_id, experiment = experiment_id, ensemble = member_id) %>%
  filter(model == unique(target_data$model) )  %>%
  select(model, experiment, ensemble, table_id, tas_zstore, psl_zstore, pr_zstore) %>%
  # cut out any rows that don't have paths to netcdfs for all variables:
  na.omit %>%
  select(model, experiment, ensemble)  %>%
  unite("acceptable", c(model, experiment, ensemble), sep="~")->
  acceptable_archive_entries_from_path_metadata

3 Generate the new ensembles

We’re actually going to generate multiple different ensembles for a target of ssp245, based on two different archives.

Archive 1 will exclude the ssp245 runs from potential matches.

Archive 2 will include the ssp245 runs - when we get to the point of targeting SSP585 runs, I think we will have to to include those runs in the archive to get at all decent generated Tgavs at the end of the century when temperature is higher. It’s still not just regurgitating the target ensembles, it will include mixing and matching from those parts.

We will also do some light exploration of the matching neighborhood size (tol argument).

# Archive data not including the target experiment:
archive_data %>% 
  # Only match to the same ESM:
  dplyr::filter(model == unique(target_data$model)) %>%  
  # Don't match to any ensemble members in the target experiment:
  dplyr::filter(experiment != unique(target_data$experiment) ) %>%
  # Final filter to make sure we only use the experiment*ensemble members with pangeo
  # path metadata for netcdfs for all varaibles of interest
  unite("acceptable", c(model, experiment, ensemble), sep="~", remove=FALSE) %>%
  filter(acceptable %in% unique(acceptable_archive_entries_from_path_metadata$acceptable)) %>%
  select(-acceptable) ->
  archive_subset1 

# what does this archive look like:
archive_subset1 %>% 
  select(experiment, ensemble) %>%
  distinct() %>% 
  group_by(experiment) %>% 
  summarize(n_complete_ensemble_members = n()) %>%
  ungroup %>%
  kable()

# Archive data including the target experiment:
archive_data %>% 
  # Only match to the same ESM:
  dplyr::filter(model == unique(target_data$model)) %>%  
  # Final filter to make sure we only use the experiment*ensemble members with pangeo
  # path metadata for netcdfs for all varaibles of interest
  unite("acceptable", c(model, experiment, ensemble), sep="~", remove=FALSE) %>%
  filter(acceptable %in% unique(acceptable_archive_entries_from_path_metadata$acceptable)) %>%
  select(-acceptable) ->
  archive_subset2 

# what does this archive look like:
archive_subset2 %>% 
  select(experiment, ensemble) %>%
  distinct() %>% 
  group_by(experiment) %>% 
  summarize(n_complete_ensemble_members = n()) %>%
  ungroup %>%
  kable()

for (draw in 1:Ndraws){
  
  # The matching to the archive without ssp245 results
  match_neighborhood(target_data = target_data,
                     archive_data = archive_subset1,
                     tol=0.1) -> matched_data
  matched_data %>%
    # Convert them to a sample of individual Tgav Recipes:
    permute_stitching_recipes(N_matches = 2000, matched_data = .,
                              archive = archive_subset1) ->
    recipes_1_0p1 
  
  # convert these recipes to something that can be ingested by python 
  # to layer in gridded netcdfs and save
  gridded_recipes_1_0p1 <- generate_gridded_recipe(recipes_1_0p1)
  write.csv(gridded_recipes_1_0p1,
            paste0('generated_ensembles/', esm_name, '/gridded_recipes_for_python_', esm_experiment,'_draw', draw, '_experiment_not_in_matching.csv'),
            row.names = FALSE)
  
  # Proceed to stitch the actual Tgavs for analysis here.
  recipes_1_0p1 %>%
    # And stitch the recipes into tgavs:
    stitch_global_mean(recipes=., data = tgav_data) %>%
    # add some identifying info
    mutate(tol = 0.1,
           archive = paste0('without_', esm_experiment)) ->
    rslts_1_0p1

  # Check for collpase in these results
    rslts_1_0p1 %>%
    select(year, value, stitching_id) %>%
    spread(year, value) ->
    wide

    # Collapse would have occured if any of these trajectories go through the same value
    # in the same year -> if the length of distinct values in any column is less than
    # the number of stitching_id's
    counts <- apply(wide, 2, length)

    print(paste('Archive without', esm_experiment,'; tol = 0.1: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))

    

  # The matching to the archive with target experiment included results
  match_neighborhood(target_data = target_data,
                     archive_data = archive_subset2,
                     tol=0.1) %>%
    # Convert them to a sample of individual Tgav Recipes:
    permute_stitching_recipes(N_matches = 2000, matched_data = .,
                              archive = archive_subset2) ->
    recipes_2_0p1
  
  # convert these recipes to something that can be ingested by python 
  # to layer in gridded netcdfs and save
  gridded_recipes_2_0p1 <- generate_gridded_recipe(recipes_2_0p1)
  write.csv(gridded_recipes_2_0p1,
            paste0('generated_ensembles/', esm_name, '/gridded_recipes_for_python_', esm_experiment,'_draw', draw, '_experiment_in_matching.csv'),
            row.names = FALSE)
  
  
  recipes_2_0p1 %>%
    # And stitch the recipes into tgavs:
    stitch_global_mean(recipes=., data = tgav_data) %>%
    # add some identifying info
    mutate(tol = 0.1,
           archive = paste0('with_', esm_experiment)) ->
    rslts_2_0p1

    # Check for collpase in these results
    rslts_2_0p1 %>%
    select(year, value, stitching_id) %>%
    spread(year, value) ->
    wide

    # Collapse would have occured if any of these trajectories go through the same value
    # in the same year -> if the length of distinct values in any column is less than
    # the number of stitching_id's
    counts <- apply(wide, 2, length)

    print(paste('Archive with', esm_experiment,'; tol = 0.1: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))

    
    
    
  # save it all off  
  write.csv(bind_rows(rslts_1_0p1, 
                      rslts_2_0p1) %>%
              mutate(draw = draw),
            paste0('generated_ensembles/', esm_name, '/generated_ensembles_', esm_experiment,'_draw', draw, '.csv'),
            row.names = FALSE)
}

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

## You have requested more recipes than possible for at least one target trajectories, returning what can

## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

4 Validate - ensembles

Comparing the generated ensembles to the training ensemble with error metrics along the lines of Tebaldi et al 2020.

That is, we need the means and variances of both the target ensemble and the generated ensemble. Specifically, mean across time and ensemble members, variance across time and ensemble members.

# values of the target data
mean_target <- mean(original_data$value)
mean_target_hist <- mean(filter(original_data, year <= 2014)$value)
mean_target_future <- mean(filter(original_data, year > 2014)$value)
  
var_target <- sd(original_data$value)
var_target_hist <- sd(filter(original_data, year <= 2014)$value)
var_target_future <- sd(filter(original_data, year > 2014)$value)
  

# generated ensemble files
# read in all generated ensembles.
files <- list.files(path = paste0('generated_ensembles/', esm_name), 
                    pattern = paste0('generated_ensembles_', esm_experiment), 
                    full.names = TRUE) 

for (i in 1:length(files)){

    # Metrics for this draw #################################################################
  generated_ensemble_draw <- read.csv(files[i], stringsAsFactors = FALSE)
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'hist+future',
           mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  generated_ensemble_draw %>% 
    filter(year <= 2014) %>%
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'hist',
           mean_target = mean_target_hist,
           var_target = var_target_hist,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics_hist
  
  
  generated_ensemble_draw%>% 
    filter(year > 2014) %>%
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'future',
           mean_target = mean_target_future,
           var_target = var_target_future,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics_future
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup ->
    tmp
  
  bind_rows(tmp %>%
              left_join(metrics, by = c('tol', 'archive', 'draw')),
            tmp %>%
              left_join(metrics_hist, by = c('tol', 'archive', 'draw')) ,
            tmp %>%
              left_join(metrics_future, by = c('tol', 'archive', 'draw')))  ->
    plotting_tbl
  rm(tmp)
  
   write.csv(plotting_tbl %>%
              mutate(draw = draw),
           paste0('generated_ensembles/', esm_name ,'/metrics_', esm_experiment,'_draw', i, '.csv'), 
            row.names = FALSE)
}

# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name),
                    pattern = paste0('metrics_', esm_experiment), 
                    full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive, time_period) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)

ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  facet_wrap(~time_period, nrow =3) +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0)

5 Validate - window jumps

Pull off the years that end each window and that start each window; these are going to be years we pay particular attention to.

window_start_years <- unique(archive_subset1$start_yr)
window_end_years <- unique(archive_subset1$end_yr)

We are particularly concerned about the ‘jumps’ at each of these years being ‘wrong’. And the idea is that they would be wrong relative to both the target data and, I think, the generated ensembles in between the jump years.

So let’s work with the jumps themselves for validating. Take the year to year diff within each realization (target or generated ensemble member). We’ll look at the diff and abs(diff). We’ll visualize the jumps with a scatter view and then get into more quantitative metrics. I spent a while going back and forth between starting from fundamentals like this and tracking down a package that looks at more sophisticated methods. As I started reading through packages, I ended up deciding to stick to fundamentals because they’ll translate from R to python better and will provide a more solid understanding of what’s going on before we look at anything more sophisticated.

# read in all generated ensembles.
files <- list.files(path = paste0('generated_ensembles/', esm_name), 
                    pattern = paste0('generated_ensembles_', esm_experiment),  
                    full.names = TRUE) 

generated_tgavs <- list()
for (i in 1:length(files)){
  generated_tgavs[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_tgavs <- do.call(rbind, generated_tgavs)


# make a helper function for calculating all the jumps because you have Nyears-1 jumps and
# dply functions don't like that. 
get_jumps <- function(data){
  data.frame(left_year = data$year[1:(length(data$year)-1)],
             jump = diff(data$value),
             abs_jump = abs(diff(data$value))) %>%
    mutate(variable = unique(data$variable),
           stitching_id = unique(data$stitching_id),
            tol = unique(data$tol),
           archive = unique(data$archive),
           draw = unique(data$draw),
           jump_year = if_else(left_year %in% window_end_years, 'window_transition', 'in_window')) ->
    x
return(x)
}

# calculate the yearly jumps
generated_list <- split(generated_tgavs, 
                        list(generated_tgavs$tol,
                        generated_tgavs$archive,
                        generated_tgavs$draw, 
                        generated_tgavs$stitching_id),
                        drop=TRUE)

lapply(generated_list, get_jumps) %>%
  do.call(rbind, .) ->
  generated_jumps

target_list <- split(original_data,
                     list(original_data$model, 
                          original_data$experiment,
                          original_data$ensemble),
                     drop = TRUE)

lapply(target_list, get_jumps) %>% 
  do.call(rbind, .) ->
  target_jumps

ggplot() + 
  geom_point(target_jumps, mapping = aes(x = left_year, y = jump), alpha = 0.5) +
  geom_point(filter(generated_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year, y = jump), alpha = 0.1, color='red') +
    facet_grid(tol~archive)

ggplot() + 
  geom_point(filter(target_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year, y = jump), alpha = 0.5) +
  geom_point(filter(generated_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year-1, y = jump), alpha = 0.1, color='red' ) +
    facet_grid(archive~tol) +
  ggtitle('note the years are staggered just for visual clarity')

6 distributions of jumps

What’s the appropriate comparison though?

1. All years target (250 jumps X 25 ensemble members) vs all years generated (250 jumps X 5-50 ensemble members depending on tol+archive X Number draws under consideration)
2. All years target (250 jumps X 25 ensemble members) vs window transition years generated (28 jumps X 5-50 ensemble members depending on tol+archive X Number draws under consideration)
3. Window transition years target (28 jumps X 25 ensemble members)vs window transition years generated (28 jumps X 5-50 ensemble members depending on tol+archive X Number draws under consideration)

From the perspective of the target data, there’s nothing special about those transition years, which basically knocks out 3. I’ll plot it for completeness but I definitely don’t think it’s the right comparison

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

focus_draw <- 1

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

focus_draw <- 2

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

# values of the target data
mean_target <- mean(target_jumps$jump)

var_target <- sd(target_jumps$jump)

  
for (drawID in 1:Ndraws){
  # Metrics for this draw #################################################################
  generated_jumps %>%
    filter(draw == drawID) ->
    generated_ensemble_draw 
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(jump), 
              var_gen = sd(jump)) %>%
    ungroup %>%
    mutate(mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup  %>%
    left_join(metrics, by = c('tol', 'archive', 'draw')) ->
    plotting_tbl
  
   write.csv(plotting_tbl %>%
              mutate(draw = drawID),
            paste0('generated_ensembles/', esm_name ,'/jump_errors_allyears_', esm_experiment, '_draw', drawID, '.csv'), 
            row.names = FALSE)
}


# Plot the first draw's errors:
plotting_tbl <- read.csv(paste0('generated_ensembles/',esm_name , '/jump_errors_allyears_',esm_experiment,'_draw1.csv'),  
                         stringsAsFactors = FALSE)

kable(plotting_tbl)

ggplot(data = plotting_tbl, aes(x = E1, y = E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle('Error metrics for single Draw 1, comparing all \nyears of generated jumps to all years of target jumps')

# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name),
                    pattern = paste0('jump_errors_allyears_', esm_experiment), 
                    full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)

ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle(paste('Aggregate errors (jumps for all years), aggregated across Ndraws ', unique(aggregate_metrics$total_draws) ))

# values of the target data
mean_target <- mean(target_jumps$jump)

var_target <- sd(target_jumps$jump)

  
for (drawID in 1:Ndraws){
  # Metrics for this draw #################################################################
  generated_jumps %>%
    filter(draw == drawID,
           jump_year=='window_transition') ->
    generated_ensemble_draw 
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(jump), 
              var_gen = sd(jump)) %>%
    ungroup %>%
    mutate(mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup  %>%
    left_join(metrics, by = c('tol', 'archive', 'draw')) ->
    plotting_tbl
  
   write.csv(plotting_tbl %>%
              mutate(draw = drawID),
            paste0('generated_ensembles/', esm_name ,'/jump_errors_transitionyears_',esm_experiment,'_draw', drawID, '.csv'), 
            row.names = FALSE)
}


# Plot the first draw's errors:
plotting_tbl <- read.csv(paste0('generated_ensembles/', esm_name ,'/jump_errors_transitionyears_',esm_experiment,'_draw1.csv'),  
                         stringsAsFactors = FALSE)

kable(plotting_tbl)

ggplot(data = plotting_tbl, aes(x = E1, y = E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle('Error metrics for single Draw 1, comparing window transition years of \ngenerated jumps to all years of target jumps')

# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name), 
                    pattern = paste0('jump_errors_transitionyears_',esm_experiment ), 
                    full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)

ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle(paste('Aggregate errors comparing window transition years of \ngenerated jumps to all years of target jumps \naggregated across Ndraws ', unique(aggregate_metrics$total_draws) ))